#include "CCell.h"
#include <QString>
#include <QMessageBox>
#include <QtDebug>


#define pi 3.14159265

int CCell::load (string fred, string fgreen, string fblue)
{
        qDebug("Started File Loading %s" ,fred.c_str());
        ifstream ifred, ifgreen, ifblue;
        ifred.open (fred.c_str (), ifstream::in);
        ifgreen.open (fgreen.c_str (), ifstream::in);
        ifblue.open (fblue.c_str (), ifstream::in);

        stringstream ssred, ssgreen, ssblue;
        string rred, ggreen, bblue;
        double r, g, b, gr;
        int j=-1;
        while (getline(ifred, rred), getline(ifgreen, ggreen), getline(ifblue, bblue))
        {
                ssred.str (rred);
                ssgreen.str (ggreen);
                ssblue.str (bblue);

                red.resize(red.size()+1);
                green.resize(green.size()+1);
                blue.resize(blue.size()+1);
                grey.resize(grey.size()+1);
                j++;

                while (ssred >> r, ssgreen >> g, ssblue >> b)
                {
                        gr=0.30*(double)r + 0.59*(double)g + 0.11*(double)b;
                        red[j].push_back(r);
                        green[j].push_back(g);
                        blue[j].push_back(b);
                        grey[j].push_back(gr);
                }
                ssred.clear ();
                ssgreen.clear ();
                ssblue.clear ();

        }

        ifred.close();
        ifgreen.close();
        ifblue.close();
        qDebug("Finished File Loading");
        return 0;
}

int CCell::genMask ()
{
    qDebug() << "Started mask generation";
    double min, max;
    min=max=grey[1][1];
    mask.resize(grey.size());
    for (int i=0; i<grey.size(); i++)
    {
            mask[i].resize(grey[i].size());
            for (int j=0; j<grey[i].size(); j++)
            {
                    if (grey[i][j]>max) max=grey[i][j];
                    if ((grey[i][j]<min)&&(grey[i][j]!=0)) min=grey[i][j];
            }
    }

    double dif=min+(max-min)*60/100;
    double razn = (max-min)/100;

    stat.resize(101);
    for (int i=0; i<grey.size(); i++)
    {
            for (int j=0; j<grey[i].size(); j++)
            {
                    if ((grey[i][j]<dif)&&(grey[i][j]!=0)) mask[i][j]=1;
                    else mask[i][j]=0;
                    int g=(grey[i][j]-min)/razn;
                    if ((g<100)&&(g>=0))
                            stat[g]++;

            }
    }
    int g=2;
    while (g--)
    for (int i=1; i<mask.size()-1; i++)
    {
        for (int j=1; j<mask[i].size()-1; j++)
        {

            if (mask[i][j])
            {
                int c=0;
                if (mask[i-1][j+1]) c++;
                if (mask[i][j+1]) c++;
                if (mask[i+1][j+1]) c++;
                if (mask[i+1][j]) c++;
                if (mask[i-1][j]) c++;
                if (mask[i+1][j-1]) c++;
                if (mask[i][j-1]) c++;
                if (mask[i-1][j-1]) c++;
                if (c<4) mask[i][j]=0;
            }
            if (!mask[i][j])
            {
                int c=0;
                if (mask[i-1][j+1]) c++;
                if (mask[i][j+1]) c++;
                if (mask[i+1][j+1]) c++;
                if (mask[i+1][j]) c++;
                if (mask[i-1][j]) c++;
                if (mask[i+1][j-1]) c++;
                if (mask[i][j-1]) c++;
                if (mask[i-1][j-1]) c++;
                if (c>6) mask[i][j]=1;
            }
        }
    }

    vector <vector <int> > mask8;
    mask8.resize(mask.size());
    for (int i=0; i< mask.size(); i++)
    {
        mask8[i].resize(mask[i].size());
    }
    for (int i=1; i< mask.size()-1; i++)
    {
        for (int j=1; j< mask[i].size()-1; j++)
        {
            int c=0;
            if (mask[i-1][j-1]) c++;
            if (mask[i-1][j]) c++;
            if (mask[i-1][j+1]) c++;
            if (mask[i][j-1]) c++;
            if (mask[i][j+1]) c++;
            if (mask[i+1][j-1]) c++;
            if (mask[i+1][j]) c++;
            if (mask[i+1][j+1]) c++;
            mask8[i][j]=c;
        }
    }

    for (int i=1; i< mask.size()-1; i++)
    {
        for (int j=1; j< mask[i].size()-1; j++)
        {
            if (mask8[i][j]!=8)
                mask[i][j]=0;
        }
    }



    vector <int> sizes;
    mask2.resize(mask.size());
    for (int i=0; i<mask.size(); i++)
    {
            mask2[i].resize(mask[i].size());
    }

    int c1=1;
    for (int i=0; i<mask.size(); i++)
    {
        for (int j=0; j<mask[i].size(); j++)
        {
            if ((mask[i][j])&&(mask2[i][j]==0))
            {
                int t=recursize (i, j, c1);
                sizes.push_back(t);
                c1++;
            }
        }
    }
    int maxv=0, maxc=0;
    for (int i=0; i<sizes.size(); i++)
    {
        if (sizes[i]>maxv)
        {
            maxv=sizes[i];
            maxc=i;
        }
    }
    maxc++;

    for (int i=0; i<mask.size(); i++)
        for (int j=0; j<mask[i].size(); j++)
        {
            if (mask2[i][j]==maxc)
                mask[i][j]=1;
            else mask[i][j]=0;
        }
    qDebug() << "Finished mask generation";
    return 0;
}

int CCell::recursize (int x, int y, int c1)
{
    if ((mask2[x][y]!=0))
        return 0;

    int res=0;
    if (mask[x][y]==1)
    {
        mask2[x][y]=c1;
        res = 1;
        if(x-1>=0)
            res+=recursize(x-1, y, c1);
        if(x+1<mask.size())
            res+=recursize(x+1, y, c1);
        if(y+1<mask[x].size())
            res+=recursize (x, y+1, c1);
        if(y-1>=0)
            res+=recursize(x, y-1, c1);
        /*if ((y-1>=0)&&(x-1>=0))
            res+=recursize(x, y-1, c1);
        if ((y-1>=0)&&(x+1<mask.size()))
            res+=recursize(x, y-1, c1);
        if ((y+1<mask[x].size())&&(x-1>=0))
            res+=recursize (x, y+1, c1);
        if ((y+1<mask[x].size())&&(x+1<mask.size()))
            res+=recursize (x, y+1, c1);
        */
    }
    return res;
}

int CCell::load (string fgrey)
{
        ifstream ifgrey;
        ifgrey.open (fgrey.c_str (), ifstream::in);

        stringstream ssgrey;
        string ggrey;
        double gr;
        int j=-1;
        while (getline(ifgrey, ggrey))
        {
                ssgrey.str (ggrey);
                grey.resize(grey.size()+1);
                j++;

                while (ssgrey >> gr)
                {
                        grey[j].push_back(gr);
                }
                ssgrey.clear ();
        }
        ifgrey.close();
        return 0;
}

int CCell::area ()
{
        int A=0;
        for (int i=0; i<mask.size(); i++)
                for (int j=0; j<mask[i].size(); j++)
                        A += mask[i][j];
        return A;
}

double CCell::x_centr ()
{
        double x_c=0;
        int ar = area();
        for (int i=0;i<mask.size(); i++)
                for (int j=0; j<mask[i].size(); j++)
                        x_c += (i*mask[i][j]*1.0)/ar;
        return x_c;
}

double CCell::y_centr ()
{
        double y_c=0;
        int ar=area();
        for (int i=0;i<mask.size(); i++)
                for (int j=0; j<mask[i].size(); j++)
                        y_c += (j*mask[i][j]*1.0)/ar;
        return y_c;
}

/*int CCell::N ()
{
        return grey.size()*grey[0].size();
}*/

double CCell::mean_radius ()
{
        N=0;
        double r=0;
        bool k;
        double x_c, y_c;
        x_c=x_centr();
        y_c=y_centr();
        for (int i=1;i<mask.size()-1; i++)
        {
                k=false;
                for (int j=1; j<mask[i].size()-1; j++)
                {
                        if ((mask[i][j]==1)&&(!k))
                                if ((mask[i][j-1]==1)&&(mask[i][j+1]==1)&&(mask[i-1][j]==1)&&(mask[i+1][j]==1)&&
                                        (mask[i-1][j-1]==1)&&(mask[i+1][j-1]==1)&&(mask[i-1][j+1]==1)&&(mask[i+1][j+1]==1))
                                {
                                        r += sqrt((i-x_c)*(i-x_c)+(j-y_c)*(j-y_c));
                                        N++;
                                        k=true;
                                }
                        if ((mask[i][j]==0)&&(k))
                        {
                                j--;
                                while (mask[i][j])
                                    if ((mask[i][j-1]==1)&&(mask[i][j+1]==1)&&(mask[i-1][j]==1)&&(mask[i+1][j]==1)&&
                                            (mask[i-1][j-1]==1)&&(mask[i+1][j-1]==1)&&(mask[i-1][j+1]==1)&&(mask[i+1][j+1]==1))
                                    {
                                            r += sqrt((i-x_c)*(i-x_c)+(j-y_c)*(j-y_c));
                                            N++;
                                            k=false;
                                            break;
                                    }
                                    else j--;
                                if (!k) break;

                        }
                }
        }
        return r/N;
}

double CCell::max_radius ()
{
        double r=0, max_r=0;
        bool k;
        double x_c, y_c;
        x_c=x_centr();
        y_c=y_centr();
        for (int i=1;i<mask.size()-1; i++)
        {
                k=false;
                for (int j=1; j<mask[i].size()-1; j++)
                {
                        if ((mask[i][j]==1)&&(!k))
                                if ((mask[i][j-1]==1)&&(mask[i][j+1]==1)&&(mask[i-1][j]==1)&&(mask[i+1][j]==1)&&
                                        (mask[i-1][j-1]==1)&&(mask[i+1][j-1]==1)&&(mask[i-1][j+1]==1)&&(mask[i+1][j+1]==1))
                                {
                                        r =sqrt((i-x_c)*(i-x_c)+(j-y_c)*(j-y_c));
                                        if (max_r<r)
                                                max_r=r;
                                        k=true;
                                }
                        if ((mask[i][j]==0)&&(k))
                        {
                                j--;
                                while (mask[i][j])
                                if ((mask[i][j-1]==1)&&(mask[i][j+1]==1)&&(mask[i-1][j]==1)&&(mask[i+1][j]==1)&&
                                        (mask[i-1][j-1]==1)&&(mask[i+1][j-1]==1)&&(mask[i-1][j+1]==1)&&(mask[i+1][j+1]==1))
                                {
                                        r =sqrt((i-x_c)*(i-x_c)+(j-y_c)*(j-y_c));
                                        if (max_r<r)
                                                max_r=r;
                                        k=false;
                                        break;
                                        }
                                        else j--;
                                        if (!k) break;

                        }
                }
        }
        return max_r;
}

double CCell::min_radius ()
{
        double r=0, min_r=2*mask.size();
        bool k;
        double x_c, y_c;
        x_c=x_centr();
        y_c=y_centr();
        for (int i=1;i<mask.size()-1; i++)
        {
                k=false;
                for (int j=1; j<mask[i].size()-1; j++)
                {
                        if ((mask[i][j]==1)&&(!k))
                                if ((mask[i][j-1]==1)&&(mask[i][j+1]==1)&&(mask[i-1][j]==1)&&(mask[i+1][j]==1)&&
                                        (mask[i-1][j-1]==1)&&(mask[i+1][j-1]==1)&&(mask[i-1][j+1]==1)&&(mask[i+1][j+1]==1))
                                {
                                        r =  sqrt((i-x_c)*(i-x_c)+(j-y_c)*(j-y_c));
                                        if (min_r>r)
                                                min_r=r;
                                        k=true;
                                }
                        if ((mask[i][j]==0)&&(k))
                        {
                                j--;
                                while (mask[i][j])
                                if ((mask[i][j-1]==1)&&(mask[i][j+1]==1)&&(mask[i-1][j]==1)&&(mask[i+1][j]==1)&&
                                        (mask[i-1][j-1]==1)&&(mask[i+1][j-1]==1)&&(mask[i-1][j+1]==1)&&(mask[i+1][j+1]==1))
                                {
                                    r = sqrt((i-x_c)*(i-x_c)+(j-y_c)*(j-y_c));
                                    if (min_r>r)
                                            min_r=r;
                                    k=false;
                                    break;
                                }
                                else j--;
                                if (!k) break;

                        }
                }
        }
        return min_r;
}

double CCell::var_radius ()
{
        N=0;
        double r, var_r=0, mean_r=mean_radius();
        bool k;
        double x_c, y_c;
        x_c=x_centr();
        y_c=y_centr();
        for (int i=1;i<mask.size()-1; i++)
        {
                k=false;
                for (int j=1; j<mask[i].size()-1; j++)
                {
                        if ((mask[i][j]==1)&&(!k))
                                if ((mask[i][j-1]==1)&&(mask[i][j+1]==1)&&(mask[i-1][j]==1)&&(mask[i+1][j]==1)&&
                                        (mask[i-1][j-1]==1)&&(mask[i+1][j-1]==1)&&(mask[i-1][j+1]==1)&&(mask[i+1][j+1]==1))
                                {
                                        r = sqrt((i-x_c)*(i-x_c)+(j-y_c)*(j-y_c));
                                        var_r+= pow(r-mean_r,2);
                                        N++;
                                        k=true;
                                }
                        if ((mask[i][j]==0)&&(k))
                        {
                                j--;
                                while (mask[i][j])
                                if ((mask[i][j-1]==1)&&(mask[i][j+1]==1)&&(mask[i-1][j]==1)&&(mask[i+1][j]==1)&&
                                        (mask[i-1][j-1]==1)&&(mask[i+1][j-1]==1)&&(mask[i-1][j+1]==1)&&(mask[i+1][j+1]==1))
                                {
                                        r += r + sqrt((i-x_c)*(i-x_c)+(j-y_c)*(j-y_c));
                                        var_r+= pow(r-mean_r,2);
                                        N++;
                                        k=false;
                                        break;
                                        }
                                        else j--;
                                        if (!k) break;

                        }
                }
        }
        return var_r/(N-1);
}

double CCell::sphericity  ()
{
        return min_radius()/max_radius();
}

double CCell::Xmoment2 ()
{
        double sum1=0;
        double XM2=0;
        for (int i=1; i<mask.size()-1;i++)
                for (int j=1; j<mask[i].size()-1; j++)
                {
                        sum1=0;
                        for (int k=1; k<mask.size()-1;k++)
                                sum1+=((double)mask[i][j]*k)/(mask.size()-1);
                        XM2+=i-sum1;
                }
        return XM2;
}

double CCell::Ymoment2 ()
{
        double sum2=0;
        double YM2=0;
        for (int i=1; i<mask.size()-1;i++)
                for (int j=1; j<mask[i].size()-1; j++)
                {
                        sum2=0;
                        for (int k=1; k<mask[i].size()-1;k++)
                                sum2+=((double)mask[i][j]*k)/(mask[i].size()-1);
                        YM2+=pow((j-sum2),2);
                }

        return YM2;
}

double CCell::XYcrossmoment2 ()
{
        double sum1=0, sum2=0;
        double XYcM2=0;
        for (int i=1; i<mask.size()-1;i++)
                for (int j=1; j<mask[i].size()-1; j++)
                {
                        sum1=sum2=0;
                        for (int k=1; k<mask.size()-1;k++)
                                sum1+=((double)mask[i][j]*k)/(mask.size()-1);
                        for (int q=1; q<mask[i].size()-1;q++)
                                sum2+=((double)mask[i][j]*q)/(mask[i].size()-1);
                        XYcM2+=(i-sum1)*(j-sum2);
                }

        return XYcM2;
}

double CCell::eccentricity ()
{
        double D=0, A1=0, A2=0, ecc=0;
        D=sqrt(pow((Xmoment2()+Ymoment2()),2)-4*(pow((XYcrossmoment2()),2)- Xmoment2()*Ymoment2() ));
        /*A1=(Xmoment2()+Ymoment2()-D)/2;
        A2=(Xmoment2()+Ymoment2()+D)/2;
        ecc=sqrt(A1/A2);*/

        return D;
}

double CCell::inertia_shape ()
{
        double ine_s=0;
        double x_c, y_c, ar;
        x_c = x_centr();
        y_c = y_centr();
        ar = area();
        for (int i=1; i<mask.size()-1; i++)
                for (int j=1; j<mask[i].size()-1; j++)
                        ine_s+=2*pi*(pow((i-x_c),2)+pow((j-y_c),2))/(pow(((double)ar),2));

        return ine_s;
}

double CCell::compactness ()
{
        double comp=0;
        double ar=area();
        for (int i=1; i<mask.size()-1; i++)
                comp=((mask.size()-1)+(mask[i].size()-1)-2)*2/(4*pi*ar);

        return comp;
}

double CCell::cell_orient ()
{
        double cell_or=0;
        double D=0, A1=0, A2=0;
        D=sqrt(pow((Xmoment2()+Ymoment2()),2)-4*(pow((XYcrossmoment2()),2)- Xmoment2()*Ymoment2() ));
        A1=(Xmoment2()+Ymoment2()-D)/2;
        A2=(Xmoment2()+Ymoment2()+D)/2;
        cell_or=(180/pi)*((pi/2)+atan((A1-Ymoment2())/XYcrossmoment2()));

        return cell_or;
}

double CCell::Ib ()
{
    double res =0;
    int count=0;
    for (int i=0; i<mask.size(); i++)
    {
        for (int j=0; j<mask[i].size(); j++)
        {
            if (mask[i][j]==0)
            {
                res+=grey[i][j];
                count++;
            }
        }
    }

    return res/count;
}



double CCell::DNA_Index ()
 {
    double DNA_Am=0;
    double backIntensity = Ib();
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            DNA_Am+=(log(backIntensity)-log(grey[i][j]))*mask[i][j];
    /*QMessageBox msgBox;
    msgBox.setText(QString::number(DNA_Am, 'f', 4));
    msgBox.setInformativeText(QString::number(DNA_Am, 'f', 4));
    msgBox.setStandardButtons(QMessageBox::Ok);
    msgBox.setDefaultButton(QMessageBox::Ok);
    int ret = msgBox.exec();*/
    return DNA_Am;
 }

double CCell::OD_max ()
{
    double OD_m=0;
    double max=-255;
    double backIntensity = Ib();
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            max = max>(log(backIntensity)-log(grey[i][j])) ? max:(log(backIntensity)-log(grey[i][j]));
    return max*100;
}

double CCell::OD_var ()
{
    double OD_v=0;
    qDebug("\tbackIntensity.clacs");
    double backIntensity = Ib();
    qDebug("\tDNA_Index.clacs");
    double DNA_Am = DNA_Index ();
    qDebug("\tarea.clacs");
    double A = area();
    qDebug("\t\tfinal.clacs");
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            OD_v+=(pow(((log(backIntensity)-log(grey[i][j]))*mask[i][j]-(DNA_Am/A)),2)/((A-1)*pow((DNA_Am/A),2)));

    /*QMessageBox msgBox;
    msgBox.setText(QString::number(OD_v, 'f', 4));
    //msgBox.setInformativeText("Do you want to save your changes?");
    msgBox.setStandardButtons(QMessageBox::Ok);
    msgBox.setDefaultButton(QMessageBox::Ok);
    int ret = msgBox.exec();*/
    return OD_v;
}

double CCell::OD_skew()
{
    double OD_sk=0;
    double backIntensity = Ib();
    double DN_IN = DNA_Index ()/area ();
    double a=area();
    double OD_s=0;
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            OD_s+=pow(((log(backIntensity)-log(grey[i][j]))*mask[i][j]-(DN_IN)),2);
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            OD_sk+=(pow(((log(backIntensity)-log(grey[i][j]))*mask[i][j]-(DN_IN)),3)/((a-1)*pow(OD_s,(double)(3/2))));
    return OD_sk;
}

double CCell::lowDNAarea ()
{
    double backIntensity = Ib();
    int Alow=0;
    for (int i=0; i<mask.size(); i++)
            for (int j=0; j<mask[i].size(); j++)
                 Alow += ((log(backIntensity)-log(grey[i][j]))>36 ? mask[i][j]:0);

    return Alow/area();
}

double CCell::hiDNAarea ()
{
    double backIntensity = Ib();
    int Ahi=0;
    for (int i=0; i<mask.size(); i++)
            for (int j=0; j<mask[i].size(); j++)
                 Ahi += ((log(backIntensity)-log(grey[i][j]))<18 ? mask[i][j]:0);

    return Ahi/area();
}

double CCell::lowDNAamnt ()
{
    double lDNAa=0;
    double backIntensity = Ib();
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            lDNAa+=(log(backIntensity)-log(grey[i][j]))>36 ? ((log(backIntensity)-log(grey[i][j]))*mask[i][j]):0;
    return lDNAa/DNA_Index();
}

double CCell::hiDNAamnt ()
{
    double hDNAa=0;
    double backIntensity = Ib();
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            hDNAa+=(log(backIntensity)-log(grey[i][j]))<18 ? ((log(backIntensity)-log(grey[i][j]))*mask[i][j]):0;
    return hDNAa/DNA_Index();
}

double CCell::hi_av_dst()
{
    double ine_s=0;
    double x_c, y_c, ar;
    double backIntensity = Ib();
    int Ahi=0;
    for (int i=0; i<mask.size(); i++)
            for (int j=0; j<mask[i].size(); j++)
                 Ahi += ((log(backIntensity)-log(grey[i][j]))<18 ? mask[i][j]:0);
    x_c = x_centr();
    y_c = y_centr();
    ar = area();
    double m_r=mean_radius();
    for (int i=1; i<mask.size()-1; i++)
            for (int j=1; j<mask[i].size()-1; j++)
                ine_s+=log(backIntensity)-log(grey[i][j])<18 ? (pow((pow((i-x_c),2)+pow((j-y_c),2)),(double)(1/2))*mask[i][j])/(Ahi*m_r):0;

    return ine_s;
}

double CCell::mh_av_dst()
{
    double ine_s=0;
    double x_c, y_c, ar;
    double backIntensity = Ib();
    int Ahi=0;
    for (int i=0; i<mask.size(); i++)
            for (int j=0; j<mask[i].size(); j++)
                 Ahi += ((log(backIntensity)-log(grey[i][j]))<18 ? mask[i][j]:0);
    int Amed=0;
    for (int i=0; i<mask.size(); i++)
            for (int j=0; j<mask[i].size(); j++)
                 Amed += ((((log(backIntensity)-log(grey[i][j]))>18)&& ((log(backIntensity)-log(grey[i][j]))<36)) ? mask[i][j]:0);
    x_c = x_centr();
    y_c = y_centr();
    ar = area();
    double m_r=mean_radius();
    double OD1=0;
    double OD2=0;
    for (int i=1; i<mask.size()-1; i++)
            for (int j=1; j<mask[i].size()-1; j++)
                OD1+=log(backIntensity)-log(grey[i][j])<18 ? (pow((pow((i-x_c),2)+pow((j-y_c),2)),(double)(1/2))*mask[i][j]):0;
    for (int i=1; i<mask.size()-1; i++)
            for (int j=1; j<mask[i].size()-1; j++)
                OD2+=log(backIntensity)-log(grey[i][j])>36 ? (pow((pow((i-x_c),2)+pow((j-y_c),2)),(double)(1/2))*mask[i][j]):0;
    ine_s=(OD1+OD2)/((Amed+Ahi)*m_r);
    return ine_s;

}
void  CCell::getDelta()
{
    delta.resize(grey.size());
    for (int i=0; i<delta.size(); i++)
        delta[i].resize(grey[i].size());
    for (int i=1; i<grey.size()-1; i++)
    {
        for (int j=1; j<grey[i].size()-1; j++)
        {
            double td=0;
            td+=fabs(grey[i][j]-grey[i-1][j-1]);
            td+=fabs(grey[i][j]-grey[i-1][j]);
            td+=fabs(grey[i][j]-grey[i-1][j+1]);
            td+=fabs(grey[i][j]-grey[i][j-1]);
            td+=fabs(grey[i][j]-grey[i][j+1]);
            td+=fabs(grey[i][j]-grey[i+1][j-1]);
            td+=fabs(grey[i][j]-grey[i+1][j]);
            td+=fabs(grey[i][j]-grey[i+1][j+1]);
            td/=8;
            delta[i][j] = td;
        }
    }
}

double CCell::entropy()
{
    double ent=0;
    for (int i=1; i<grey.size()-1; i++)
        for (int j=1; j<grey[i].size()-1; j++)
            ent+=delta[i][j]*log(delta[i][j]);
    return ent;
}

double CCell::correlation ()
{
    double Ii=0;
    int count=0;
    for (int i=1; i<mask.size()-1; i++)
    {     for (int j=1; j<mask[i].size()-1; j++)
        {Ii+= grey[i][j]<40 ? grey[i][j]:0;
        count++;}
    }
    Ii=Ii/count;

    double cor=0;
    for (int i=1; i<grey.size()-1; i++)
        for (int j=1; j<grey[i].size()-1; j++)
            cor+=(i-Ii)*(j-Ii)*delta[i][j];
    return cor;

}

double CCell::homogeneity ()
{
    double hom=0;
    for (int i=1; i<grey.size()-1; i++)
        for (int j=1; j<grey[i].size()-1; j++)
            hom+=(1/(1+pow((i-j),2.0)))*delta[i][j];
    return hom;

}

double CCell::den_lit_spot ()
{
    double max=-255;
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            max = max>grey[i][j] ? max:grey[i][j];

    double dls=0;
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            dls+=grey[i][j]==max ? 1:0;
    return dls/area();
}

double CCell::renge_average ()
{
    double max=-255;
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            max = max>grey[i][j] ? max:grey[i][j];

    double min=-255;
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            min = min<grey[i][j] ? min:grey[i][j];

    double dmax=0;
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            dmax+=grey[i][j]==max ? 1:0;

    double dmin=0;
    for (int i=1; i<mask.size()-1; i++)
        for (int j=1; j<mask[i].size()-1; j++)
            dmin+=grey[i][j]==min ? 1:0;

    return (max/dmax-min/dmin)*100;
}

double CCell::fractal1_area ()
{
    double f1_a=0;
    double backIntensity = Ib();
    for (int i=2; i<mask.size()-1; i++)
        for (int j=2; j<mask[i].size()-1; j++)
            f1_a+=(abs((log(backIntensity)-log(grey[i][j]))-(log(backIntensity)-log(grey[i][j-1])))+abs((log(backIntensity)-log(grey[i][j]))-(log(backIntensity)-log(grey[i-1][j])))+1)*mask[i][j];
    return f1_a;
}

double CCell::fractal2_area ()
{
    double f2_a=0;
    double backIntensity = Ib();
    for (int i=2; i<(mask.size()-1)/2; i++)
        for (int j=2; j<(mask[i].size()-1)/2; j++)
            f2_a+=(abs((log(backIntensity)-log(grey[i][j]))-(log(backIntensity)-log(grey[i][j-1])))+abs((log(backIntensity)-log(grey[i][j]))-(log(backIntensity)-log(grey[i-1][j])))+1)*mask[i][j];
    return f2_a;
}

double CCell::fractal_dimen ()
{
    return (log(fractal1_area())-log(fractal2_area()))/log(2.0);
}

double CCell::grey_level()
{
    double g_l1=0;
    double g_l2=0;
     double g_l3=0;
    double x_c = x_centr();
    double y_c = y_centr();
    double count1=0;
    map <double, int> levelCount;

    for (int i=1; i<mask.size()-1; i++)
         for (int j=1; j<mask[i].size()-1; j++)
             levelCount[grey[i][j]]++;

    for (int i=1; i<levelCount.size(); i++)
    {
        {
            for (int j=1; j<8; j++)
                    g_l1+=pow((pow((i-x_c),2)+pow((j-y_c),2)),(double)(1/2));
        }
        g_l2+=pow(g_l1,2);
    }

    for (int i=1; i<levelCount.size()-1; i++)
         for (int j=1; j<8; j++)
                    g_l3+=pow((pow((i-x_c),2)+pow((j-y_c),2)),(double)(1/2));

    return g_l2/g_l3;
}

double CCell::run_length()
{
    double g_l1=0;
    double g_l2=0;
     double g_l3=0;
    double x_c = x_centr();
    double y_c = y_centr();
    double count1=0;
    map <double, int> levelCount;

    for (int i=1; i<mask.size()-1; i++)
         for (int j=1; j<mask[i].size()-1; j++)
             levelCount[grey[i][j]]++;

    for (int i=1; i<8; i++)
    {
        {
            for (int j=1; j<levelCount.size(); j++)
                    g_l1+=pow((pow((i-x_c),2)+pow((j-y_c),2)),(double)(1/2));
        }
        g_l2+=pow(g_l1,2);
    }

    for (int i=1; i<levelCount.size()-1; i++)
         for (int j=1; j<8; j++)
                    g_l3+=pow((pow((i-x_c),2)+pow((j-y_c),2)),(double)(1/2));

    return g_l2/g_l3;
}

double CCell::run_persent()
{
    double g_l2=0;
    double x_c = x_centr();
    double y_c = y_centr();
    map <double, int> levelCount;

    for (int i=1; i<mask.size()-1; i++)
         for (int j=1; j<mask[i].size()-1; j++)
             levelCount[grey[i][j]]++;


    for (int i=1; i<levelCount.size()-1; i++)
         for (int j=1; j<8; j++)
                    g_l2+=pow((pow((i-x_c),2)+pow((j-y_c),2)),(double)(1/2));

    return g_l2/area();
}

int CCell::Ozu(vector<vector<double> > grey)
{
    vector<int> image;
    for (int i=1; i<grey.size(); i++)
        for (int j=1; j<grey[i].size(); j++)
        image.push_back((int)grey[i][j]);
    int min = image[0];
    int max = image[0];
    for (int i = 1; i < image.size()-1; i++)
        {
            int value = image[i];
            if (value < min)
            min = value;
            if (value > max)
            max = value;
        }

    int histSize = max - min + 1;
    int* hist = new int[histSize];

    for (int t = 0; t < histSize; t++)
    hist[t] = 0;

    for (int i = 0; i < image.size()-1; i++)
    hist[image[i] - min]++;

    int m = 0;
    int n = 0;

    for (int t = 0; t <= max - min; t++)
    {
     m += t * hist[t];
     n += hist[t];
          }

          float maxSigma = -1;
          int threshold = 0;

          int alpha1 = 0;
          int beta1 = 0;

          for (int t = 0; t < max - min; t++)
          {
            alpha1 += t * hist[t];
            beta1 += hist[t];


            float w1 = (float)beta1 / n;

            float a = (float)alpha1 / beta1 - (float)(m - alpha1) / (n - beta1);


            float sigma = w1 * (1 - w1) * a * a;

            if (sigma > maxSigma)
            {
              maxSigma = sigma;
              threshold = t;
            }
          }


          threshold += min;


          return threshold;



}

vector<double> CCell::getSignals () {
    vector <double> res;
    res.push_back(this->area());
    qDebug("Started calculations: mean_radius");
    res.push_back(this->mean_radius());
    qDebug("Started calculations: max_radius");
    res.push_back(this->max_radius());
    qDebug("Started calculations: var_radius");
    res.push_back(this->var_radius());
    qDebug("Started calculations: sphericity()");
    res.push_back(this->sphericity());
    qDebug("Started calculations: inertia_shape()");
    res.push_back(this->inertia_shape());
    qDebug("Started calculations: compactness()");
    res.push_back(this->compactness());
    qDebug("Started calculations: OD_var()");
    res.push_back(this->OD_var());
    qDebug("Started calculations: OD_skew()");
    res.push_back(this->OD_skew());
    qDebug("Started calculations: range_avarege()");
    res.push_back(this->renge_average());
    qDebug("Started calculations: OD_max()");
    res.push_back(this->OD_max());
    qDebug("Started calculations: den_lit_spot()");
    res.push_back(this->den_lit_spot());
    qDebug("Started calculations: lowDNaarea()");
    res.push_back(this->lowDNAarea());
    qDebug("Started calculations: hiDNaarea");
    res.push_back(this->hiDNAarea());
    qDebug("Started calculations: lowDNaamnt");
    res.push_back(this->lowDNAamnt());
    qDebug("Started calculations: hiDNaamnt");
    res.push_back(this->hiDNAamnt());
    qDebug("Started calculations: hi_av_dst");
    res.push_back(this->hi_av_dst());
    qDebug("Started calculations: mh_av_dst");
    res.push_back(this->mh_av_dst());
    qDebug("Started calculations: correlation");
    res.push_back(this->correlation());
    qDebug("Started calculations: homogeneity");
    res.push_back(this->homogeneity());
    qDebug("Started calculations: entropy");
    res.push_back(this->entropy());
    qDebug("Started calculations: fractal_dimen");
    res.push_back(this->fractal_dimen());
    qDebug("Started calculations: DNa_Index");
    res.push_back(this->DNA_Index());
    qDebug("Started calculations: run_persent");
    res.push_back(this->run_persent());
    qDebug("Finished");
    return res;
}
